home *** CD-ROM | disk | FTP | other *** search
/ The CICA Windows Explosion! / The CICA Windows Explosion! - Disc 2.iso / programr / dpmigcc5.zip / RSX / SOURCE / FPU-EMU / FPU_TRIG.C < prev    next >
C/C++ Source or Header  |  1994-05-27  |  41KB  |  1,742 lines
  1. /*---------------------------------------------------------------------------+
  2.  |  fpu_trig.c                                                               |
  3.  |                                                                           |
  4.  | Implementation of the FPU "transcendental" functions.                     |
  5.  |                                                                           |
  6.  | Copyright (C) 1992,1993,1994                                              |
  7.  |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
  8.  |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
  9.  |                                                                           |
  10.  |                                                                           |
  11.  +---------------------------------------------------------------------------*/
  12.  
  13. #include "fpu_system.h"
  14. #include "exception.h"
  15. #include "fpu_emu.h"
  16. #include "status_w.h"
  17. #include "control_w.h"
  18. #include "reg_constant.h"    
  19.  
  20.  
  21. static void rem_kernel(unsigned long long st0, unsigned long long *y,
  22.                unsigned long long st1,
  23.                unsigned long long q, int n);
  24.  
  25. #define BETTER_THAN_486
  26.  
  27. #define FCOS  4
  28. #define FPTAN 1
  29.  
  30.  
  31. /* Used only by fptan, fsin, fcos, and fsincos. */
  32. /* This routine produces very accurate results, similar to
  33.    using a value of pi with more than 128 bits precision. */
  34. /* Limited measurements show no results worse than 64 bit precision
  35.    except for the results for arguments close to 2^63, where the
  36.    precision of the result sometimes degrades to about 63.9 bits */
  37. static int trig_arg(FPU_REG *X, int even)
  38. {
  39.   FPU_REG tmp;
  40.   unsigned long long q;
  41.   int old_cw = control_word, saved_status = partial_status;
  42.  
  43.   if ( X->exp >= EXP_BIAS + 63 )
  44.     {
  45.       partial_status |= SW_C2;     /* Reduction incomplete. */
  46.       return -1;
  47.     }
  48.  
  49.   control_word &= ~CW_RC;
  50.   control_word |= RC_CHOP;
  51.  
  52.   reg_div(X, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f);
  53.   round_to_int(&tmp);  /* Fortunately, this can't overflow
  54.               to 2^64 */
  55.   q = significand(&tmp);
  56.   if ( q )
  57.     {
  58.       rem_kernel(significand(X),
  59.          &significand(&tmp),
  60.          significand(&CONST_PI2),
  61.          q, X->exp - CONST_PI2.exp);
  62.       tmp.exp = CONST_PI2.exp;
  63.       normalize(&tmp);
  64.       reg_move(&tmp, X);
  65.     }
  66.  
  67.   if ( even == FPTAN )
  68.     {
  69.       if ( ((X->exp >= EXP_BIAS) ||
  70.         ((X->exp == EXP_BIAS-1)
  71.          && (X->sigh >= 0xc90fdaa2))) ^ (q & 1) )
  72.     even = FCOS;
  73.       else
  74.     even = 0;
  75.     }
  76.  
  77.   if ( (even && !(q & 1)) || (!even && (q & 1)) )
  78.     {
  79.       reg_sub(&CONST_PI2, X, X, FULL_PRECISION);
  80. #ifdef BETTER_THAN_486
  81.       /* So far, the results are exact but based upon a 64 bit
  82.      precision approximation to pi/2. The technique used
  83.      now is equivalent to using an approximation to pi/2 which
  84.      is accurate to about 128 bits. */
  85.       if ( (X->exp <= CONST_PI2extra.exp + 64) || (q > 1) )
  86.     {
  87.       /* This code gives the effect of having p/2 to better than
  88.          128 bits precision. */
  89.       significand(&tmp) = q + 1;
  90.       tmp.exp = EXP_BIAS + 63;
  91.       tmp.tag = TW_Valid;
  92.       normalize(&tmp);
  93.       reg_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION);
  94.       reg_add(X, &tmp,  X, FULL_PRECISION);
  95.       if ( X->sign == SIGN_NEG )
  96.         {
  97.           /* CONST_PI2extra is negative, so the result of the addition
  98.          can be negative. This means that the argument is actually
  99.          in a different quadrant. The correction is always < pi/2,
  100.          so it can't overflow into yet another quadrant. */
  101.           X->sign = SIGN_POS;
  102.           q++;
  103.         }
  104.     }
  105. #endif BETTER_THAN_486
  106.     }
  107. #ifdef BETTER_THAN_486
  108.   else
  109.     {
  110.       /* So far, the results are exact but based upon a 64 bit
  111.      precision approximation to pi/2. The technique used
  112.      now is equivalent to using an approximation to pi/2 which
  113.      is accurate to about 128 bits. */
  114.       if ( ((q > 0) && (X->exp <= CONST_PI2extra.exp + 64)) || (q > 1) )
  115.     {
  116.       /* This code gives the effect of having p/2 to better than
  117.          128 bits precision. */
  118.       significand(&tmp) = q;
  119.       tmp.exp = EXP_BIAS + 63;
  120.       tmp.tag = TW_Valid;
  121.       normalize(&tmp);
  122.       reg_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION);
  123.       reg_sub(X, &tmp, X, FULL_PRECISION);
  124.       if ( (X->exp == CONST_PI2.exp) &&
  125.           ((X->sigh > CONST_PI2.sigh)
  126.            || ((X->sigh == CONST_PI2.sigh)
  127.            && (X->sigl > CONST_PI2.sigl))) )
  128.         {
  129.           /* CONST_PI2extra is negative, so the result of the
  130.          subtraction can be larger than pi/2. This means
  131.          that the argument is actually in a different quadrant.
  132.          The correction is always < pi/2, so it can't overflow
  133.          into yet another quadrant. */
  134.           reg_sub(&CONST_PI, X, X, FULL_PRECISION);
  135.           q++;
  136.         }
  137.     }
  138.     }
  139. #endif BETTER_THAN_486
  140.  
  141.   control_word = old_cw;
  142.   partial_status = saved_status & ~SW_C2;     /* Reduction complete. */
  143.  
  144.   return (q & 3) | even;
  145. }
  146.  
  147.  
  148. /* Convert a long to register */
  149. void convert_l2reg(long const *arg, FPU_REG *dest)
  150. {
  151.   long num = *arg;
  152.  
  153.   if (num == 0)
  154.     { reg_move(&CONST_Z, dest); return; }
  155.  
  156.   if (num > 0)
  157.     dest->sign = SIGN_POS;
  158.   else
  159.     { num = -num; dest->sign = SIGN_NEG; }
  160.  
  161.   dest->sigh = num;
  162.   dest->sigl = 0;
  163.   dest->exp = EXP_BIAS + 31;
  164.   dest->tag = TW_Valid;
  165.   normalize(dest);
  166. }
  167.  
  168.  
  169. static void single_arg_error(void)
  170. {
  171.   switch ( FPU_st0_tag )
  172.     {
  173.     case TW_NaN:
  174.       if ( !(FPU_st0_ptr->sigh & 0x40000000) )   /* Signaling ? */
  175.     {
  176.       EXCEPTION(EX_Invalid);
  177.       if ( control_word & CW_Invalid )
  178.         FPU_st0_ptr->sigh |= 0x40000000;      /* Convert to a QNaN */
  179.     }
  180.       break;              /* return with a NaN in st(0) */
  181.     case TW_Empty:
  182.       stack_underflow();  /* Puts a QNaN in st(0) */
  183.       break;
  184. #ifdef PARANOID
  185.     default:
  186.       EXCEPTION(EX_INTERNAL|0x0112);
  187. #endif PARANOID
  188.     }
  189. }
  190.  
  191.  
  192. static void single_arg_2_error(void)
  193. {
  194.   FPU_REG *st_new_ptr;
  195.  
  196.   switch ( FPU_st0_tag )
  197.     {
  198.     case TW_NaN:
  199.       if ( !(FPU_st0_ptr->sigh & 0x40000000) )   /* Signaling ? */
  200.     {
  201.       EXCEPTION(EX_Invalid);
  202.       if ( control_word & CW_Invalid )
  203.         {
  204.           /* The masked response */
  205.           /* Convert to a QNaN */
  206.           FPU_st0_ptr->sigh |= 0x40000000;
  207.           st_new_ptr = &st(-1);
  208.           push();
  209.           reg_move(&st(1), FPU_st0_ptr);
  210.         }
  211.     }
  212.       else
  213.     {
  214.       /* A QNaN */
  215.       st_new_ptr = &st(-1);
  216.       push();
  217.       reg_move(&st(1), FPU_st0_ptr);
  218.     }
  219.       break;              /* return with a NaN in st(0) */
  220. #ifdef PARANOID
  221.     default:
  222.       EXCEPTION(EX_INTERNAL|0x0112);
  223. #endif PARANOID
  224.     }
  225. }
  226.  
  227.  
  228. /*---------------------------------------------------------------------------*/
  229.  
  230. static void f2xm1(void)
  231. {
  232.   clear_C1();
  233.   switch ( FPU_st0_tag )
  234.     {
  235.     case TW_Valid:
  236.       {
  237.     FPU_REG rv, tmp;
  238.  
  239.     if ( FPU_st0_ptr->exp >= 0 )
  240.       {
  241.         /* For an 80486 FPU, the result is undefined. */
  242.       }
  243.     else if ( FPU_st0_ptr->exp >= -64 )
  244.       {
  245.         if ( FPU_st0_ptr->sign == SIGN_POS )
  246.           {
  247.         /* poly_2xm1(x) requires 0 < x < 1. */
  248.         poly_2xm1(FPU_st0_ptr, &rv);
  249.         reg_mul(&rv, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
  250.           }
  251.         else
  252.           {
  253.         /* poly_2xm1(x) doesn't handle negative numbers yet. */
  254.         /* So we compute z=poly_2xm1(-x), and the answer is
  255.            then -z/(1+z) */
  256.         FPU_st0_ptr->sign = SIGN_POS;
  257.         poly_2xm1(FPU_st0_ptr, &rv);
  258.         reg_mul(&rv, FPU_st0_ptr, &rv, FULL_PRECISION);
  259.         reg_add(&rv, &CONST_1, &tmp, FULL_PRECISION);
  260.         reg_div(&rv, &tmp, FPU_st0_ptr, FULL_PRECISION);
  261.         FPU_st0_ptr->sign = SIGN_NEG;
  262.           }
  263.       }
  264.     else
  265.       {
  266. #ifdef DENORM_OPERAND
  267.         if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  268.           return;
  269. #endif DENORM_OPERAND
  270.         /* For very small arguments, this is accurate enough. */
  271.         reg_mul(&CONST_LN2, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
  272.       }
  273.     set_precision_flag_up();
  274.     return;
  275.       }
  276.     case TW_Zero:
  277.       return;
  278.     case TW_Infinity:
  279.       if ( FPU_st0_ptr->sign == SIGN_NEG )
  280.     {
  281.       /* -infinity gives -1 (p16-10) */
  282.       reg_move(&CONST_1, FPU_st0_ptr);
  283.       FPU_st0_ptr->sign = SIGN_NEG;
  284.     }
  285.       return;
  286.     default:
  287.       single_arg_error();
  288.     }
  289. }
  290.  
  291.  
  292. static void fptan(void)
  293. {
  294.   FPU_REG *st_new_ptr;
  295.   int q;
  296.   char arg_sign = FPU_st0_ptr->sign;
  297.  
  298.   /* Stack underflow has higher priority */
  299.   if ( FPU_st0_tag == TW_Empty )
  300.     {
  301.       stack_underflow();  /* Puts a QNaN in st(0) */
  302.       if ( control_word & CW_Invalid )
  303.     {
  304.       st_new_ptr = &st(-1);
  305.       push();
  306.       stack_underflow();  /* Puts a QNaN in the new st(0) */
  307.     }
  308.       return;
  309.     }
  310.  
  311.   if ( STACK_OVERFLOW )
  312.     { stack_overflow(); return; }
  313.  
  314.   switch ( FPU_st0_tag )
  315.     {
  316.     case TW_Valid:
  317.  
  318.       if ( FPU_st0_ptr->exp > EXP_BIAS - 40 )
  319.     {
  320.       FPU_st0_ptr->sign = SIGN_POS;
  321.       if ( (q = trig_arg(FPU_st0_ptr, FPTAN)) != -1 )
  322.         {
  323.           reg_div(FPU_st0_ptr, &CONST_PI2, FPU_st0_ptr,
  324.               FULL_PRECISION);
  325.           poly_tan(FPU_st0_ptr, FPU_st0_ptr, q & FCOS);
  326.           FPU_st0_ptr->sign = (q & 1) ^ arg_sign;
  327.         }
  328.       else
  329.         {
  330.           /* Operand is out of range */
  331.           FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
  332.           return;
  333.         }
  334.     }
  335.       else
  336.     {
  337.       /* For a small arg, the result == the argument */
  338.       /* Underflow may happen */
  339.  
  340.       if ( FPU_st0_ptr->exp <= EXP_UNDER )
  341.         {
  342. #ifdef DENORM_OPERAND
  343.           if ( denormal_operand() )
  344.         return;
  345. #endif DENORM_OPERAND
  346.           /* A denormal result has been produced.
  347.          Precision must have been lost, this is always
  348.          an underflow. */
  349.           if ( arith_underflow(FPU_st0_ptr) )
  350.         return;
  351.         }
  352.       else
  353.         set_precision_flag_up();  /* Must be up. */
  354.     }
  355.       push();
  356.       reg_move(&CONST_1, FPU_st0_ptr);
  357.       return;
  358.       break;
  359.     case TW_Infinity:
  360.       /* The 80486 treats infinity as an invalid operand */
  361.       arith_invalid(FPU_st0_ptr);
  362.       if ( control_word & CW_Invalid )
  363.     {
  364.       st_new_ptr = &st(-1);
  365.       push();
  366.       arith_invalid(FPU_st0_ptr);
  367.     }
  368.       return;
  369.     case TW_Zero:
  370.       push();
  371.       reg_move(&CONST_1, FPU_st0_ptr);
  372.       setcc(0);
  373.       break;
  374.     default:
  375.       single_arg_2_error();
  376.       break;
  377.     }
  378. }
  379.  
  380.  
  381. static void fxtract(void)
  382. {
  383.   FPU_REG *st_new_ptr;
  384.   register FPU_REG *st1_ptr = FPU_st0_ptr;  /* anticipate */
  385.  
  386.   if ( STACK_OVERFLOW )
  387.     {  stack_overflow(); return; }
  388.   clear_C1();
  389.   if ( !(FPU_st0_tag ^ TW_Valid) )
  390.     {
  391.       long e;
  392.  
  393. #ifdef DENORM_OPERAND
  394.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  395.     return;
  396. #endif DENORM_OPERAND
  397.       
  398.       push();
  399.       reg_move(st1_ptr, FPU_st0_ptr);
  400.       FPU_st0_ptr->exp = EXP_BIAS;
  401.       e = st1_ptr->exp - EXP_BIAS;
  402.       convert_l2reg(&e, st1_ptr);
  403.       return;
  404.     }
  405.   else if ( FPU_st0_tag == TW_Zero )
  406.     {
  407.       char sign = FPU_st0_ptr->sign;
  408.       if ( divide_by_zero(SIGN_NEG, FPU_st0_ptr) )
  409.     return;
  410.       push();
  411.       reg_move(&CONST_Z, FPU_st0_ptr);
  412.       FPU_st0_ptr->sign = sign;
  413.       return;
  414.     }
  415.   else if ( FPU_st0_tag == TW_Infinity )
  416.     {
  417.       char sign = FPU_st0_ptr->sign;
  418.       FPU_st0_ptr->sign = SIGN_POS;
  419.       push();
  420.       reg_move(&CONST_INF, FPU_st0_ptr);
  421.       FPU_st0_ptr->sign = sign;
  422.       return;
  423.     }
  424.   else if ( FPU_st0_tag == TW_NaN )
  425.     {
  426.       if ( real_2op_NaN(FPU_st0_ptr, FPU_st0_ptr, FPU_st0_ptr) )
  427.     return;
  428.       push();
  429.       reg_move(st1_ptr, FPU_st0_ptr);
  430.       return;
  431.     }
  432.   else if ( FPU_st0_tag == TW_Empty )
  433.     {
  434.       /* Is this the correct behaviour? */
  435.       if ( control_word & EX_Invalid )
  436.     {
  437.       stack_underflow();
  438.       push();
  439.       stack_underflow();
  440.     }
  441.       else
  442.     EXCEPTION(EX_StackUnder);
  443.     }
  444. #ifdef PARANOID
  445.   else
  446.     EXCEPTION(EX_INTERNAL | 0x119);
  447. #endif PARANOID
  448. }
  449.  
  450.  
  451. static void fdecstp(void)
  452. {
  453.   clear_C1();
  454.   top--;  /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
  455. }
  456.  
  457. static void fincstp(void)
  458. {
  459.   clear_C1();
  460.   top++;  /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
  461. }
  462.  
  463.  
  464. static void fsqrt_(void)
  465. {
  466.   clear_C1();
  467.   if ( !(FPU_st0_tag ^ TW_Valid) )
  468.     {
  469.       int expon;
  470.       
  471.       if (FPU_st0_ptr->sign == SIGN_NEG)
  472.     {
  473.       arith_invalid(FPU_st0_ptr);  /* sqrt(negative) is invalid */
  474.       return;
  475.     }
  476.  
  477. #ifdef DENORM_OPERAND
  478.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  479.     return;
  480. #endif DENORM_OPERAND
  481.  
  482.       expon = FPU_st0_ptr->exp - EXP_BIAS;
  483.       FPU_st0_ptr->exp = EXP_BIAS + (expon & 1);  /* make st(0) in  [1.0 .. 4.0) */
  484.       
  485.       wm_sqrt(FPU_st0_ptr, control_word);    /* Do the computation */
  486.       
  487.       FPU_st0_ptr->exp += expon >> 1;
  488.       FPU_st0_ptr->sign = SIGN_POS;
  489.     }
  490.   else if ( FPU_st0_tag == TW_Zero )
  491.     return;
  492.   else if ( FPU_st0_tag == TW_Infinity )
  493.     {
  494.       if ( FPU_st0_ptr->sign == SIGN_NEG )
  495.     arith_invalid(FPU_st0_ptr);  /* sqrt(-Infinity) is invalid */
  496.       return;
  497.     }
  498.   else
  499.     { single_arg_error(); return; }
  500.  
  501. }
  502.  
  503.  
  504. static void frndint_(void)
  505. {
  506.   int flags;
  507.  
  508.   if ( !(FPU_st0_tag ^ TW_Valid) )
  509.     {
  510.       if (FPU_st0_ptr->exp > EXP_BIAS+63)
  511.     return;
  512.  
  513. #ifdef DENORM_OPERAND
  514.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  515.     return;
  516. #endif DENORM_OPERAND
  517.  
  518.       /* Fortunately, this can't overflow to 2^64 */
  519.       if ( (flags = round_to_int(FPU_st0_ptr)) )
  520.     set_precision_flag(flags);
  521.  
  522.       FPU_st0_ptr->exp = EXP_BIAS + 63;
  523.       normalize(FPU_st0_ptr);
  524.       return;
  525.     }
  526.   else if ( (FPU_st0_tag == TW_Zero) || (FPU_st0_tag == TW_Infinity) )
  527.     return;
  528.   else
  529.     single_arg_error();
  530. }
  531.  
  532.  
  533. static void fsin(void)
  534. {
  535.   char arg_sign = FPU_st0_ptr->sign;
  536.  
  537.   if ( FPU_st0_tag == TW_Valid )
  538.     {
  539.       FPU_REG rv;
  540.       int q;
  541.  
  542.       if ( FPU_st0_ptr->exp > EXP_BIAS - 40 )
  543.     {
  544.       FPU_st0_ptr->sign = SIGN_POS;
  545.       if ( (q = trig_arg(FPU_st0_ptr, 0)) != -1 )
  546.         {
  547.           reg_div(FPU_st0_ptr, &CONST_PI2, FPU_st0_ptr, FULL_PRECISION);
  548.  
  549.           poly_sine(FPU_st0_ptr, &rv);
  550.  
  551.           if (q & 2)
  552.         rv.sign ^= SIGN_POS ^ SIGN_NEG;
  553.           rv.sign ^= arg_sign;
  554.           reg_move(&rv, FPU_st0_ptr);
  555.  
  556.           /* We do not really know if up or down */
  557.           set_precision_flag_up();
  558.           return;
  559.         }
  560.       else
  561.         {
  562.           /* Operand is out of range */
  563.           FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
  564.           return;
  565.         }
  566.     }
  567.       else
  568.     {
  569.       /* For a small arg, the result == the argument */
  570.       /* Underflow may happen */
  571.  
  572.       if ( FPU_st0_ptr->exp <= EXP_UNDER )
  573.         {
  574. #ifdef DENORM_OPERAND
  575.           if ( denormal_operand() )
  576.         return;
  577. #endif DENORM_OPERAND
  578.           /* A denormal result has been produced.
  579.          Precision must have been lost, this is always
  580.          an underflow. */
  581.           arith_underflow(FPU_st0_ptr);
  582.           return;
  583.         }
  584.  
  585.       set_precision_flag_up();  /* Must be up. */
  586.     }
  587.     }
  588.   else if ( FPU_st0_tag == TW_Zero )
  589.     {
  590.       setcc(0);
  591.       return;
  592.     }
  593.   else if ( FPU_st0_tag == TW_Infinity )
  594.     {
  595.       /* The 80486 treats infinity as an invalid operand */
  596.       arith_invalid(FPU_st0_ptr);
  597.       return;
  598.     }
  599.   else
  600.     single_arg_error();
  601. }
  602.  
  603.  
  604. static int f_cos(FPU_REG *arg)
  605. {
  606.   char arg_sign = arg->sign;
  607.  
  608.   if ( arg->tag == TW_Valid )
  609.     {
  610.       FPU_REG rv;
  611.       int q;
  612.  
  613.       if ( arg->exp > EXP_BIAS - 40 )
  614.     {
  615.       arg->sign = SIGN_POS;
  616.       if ( (q = trig_arg(arg, FCOS)) != -1 )
  617.         {
  618.           reg_div(arg, &CONST_PI2, arg, FULL_PRECISION);
  619.           
  620.           poly_sine(arg, &rv);
  621.  
  622.           if ((q+1) & 2)
  623.         rv.sign ^= SIGN_POS ^ SIGN_NEG;
  624.           reg_move(&rv, arg);
  625.  
  626.           /* We do not really know if up or down */
  627.           set_precision_flag_down();
  628.       
  629.           return 0;
  630.         }
  631.       else
  632.         {
  633.           /* Operand is out of range */
  634.           arg->sign = arg_sign;         /* restore st(0) */
  635.           return 1;
  636.         }
  637.     }
  638.       else
  639.     {
  640. #ifdef DENORM_OPERAND
  641.       if ( (arg->exp <= EXP_UNDER) && (denormal_operand()) )
  642.         return 1;
  643. #endif DENORM_OPERAND
  644.  
  645.       setcc(0);
  646.       reg_move(&CONST_1, arg);
  647. #ifdef PECULIAR_486
  648.       set_precision_flag_down();  /* 80486 appears to do this. */
  649. #else
  650.       set_precision_flag_up();  /* Must be up. */
  651. #endif PECULIAR_486
  652.       return 0;
  653.     }
  654.     }
  655.   else if ( arg->tag == TW_Zero )
  656.     {
  657.       reg_move(&CONST_1, arg);
  658.       setcc(0);
  659.       return 0;
  660.     }
  661.   else if ( FPU_st0_tag == TW_Infinity )
  662.     {
  663.       /* The 80486 treats infinity as an invalid operand */
  664.       arith_invalid(FPU_st0_ptr);
  665.       return 1;
  666.     }
  667.   else
  668.     {
  669.       single_arg_error();  /* requires arg == &st(0) */
  670.       return 1;
  671.     }
  672. }
  673.  
  674.  
  675. static void fcos(void)
  676. {
  677.   f_cos(FPU_st0_ptr);
  678. }
  679.  
  680.  
  681. static void fsincos(void)
  682. {
  683.   FPU_REG *st_new_ptr;
  684.   FPU_REG arg;
  685.  
  686.   /* Stack underflow has higher priority */
  687.   if ( FPU_st0_tag == TW_Empty )
  688.     {
  689.       stack_underflow();  /* Puts a QNaN in st(0) */
  690.       if ( control_word & CW_Invalid )
  691.     {
  692.       st_new_ptr = &st(-1);
  693.       push();
  694.       stack_underflow();  /* Puts a QNaN in the new st(0) */
  695.     }
  696.       return;
  697.     }
  698.  
  699.   if ( STACK_OVERFLOW )
  700.     { stack_overflow(); return; }
  701.  
  702.   if ( FPU_st0_tag == TW_NaN )
  703.     {
  704.       single_arg_2_error();
  705.       return;
  706.     }
  707.   else if ( FPU_st0_tag == TW_Infinity )
  708.     {
  709.       /* The 80486 treats infinity as an invalid operand */
  710.       if ( !arith_invalid(FPU_st0_ptr) )
  711.     {
  712.       /* unmasked response */
  713.       push();
  714.       arith_invalid(FPU_st0_ptr);
  715.     }
  716.       return;
  717.     }
  718.  
  719.   reg_move(FPU_st0_ptr,&arg);
  720.   if ( !f_cos(&arg) )
  721.     {
  722.       fsin();
  723.       push();
  724.       reg_move(&arg,FPU_st0_ptr);
  725.     }
  726.  
  727. }
  728.  
  729.  
  730. /*---------------------------------------------------------------------------*/
  731. /* The following all require two arguments: st(0) and st(1) */
  732.  
  733. /* A lean, mean kernel for the fprem instructions. This relies upon
  734.    the division and rounding to an integer in do_fprem giving an
  735.    exact result. Because of this, rem_kernel() needs to deal only with
  736.    the least significant 64 bits, the more significant bits of the
  737.    result must be zero.
  738.  */
  739. static void rem_kernel(unsigned long long st0, unsigned long long *y,
  740.                unsigned long long st1,
  741.                unsigned long long q, int n)
  742. {
  743.   unsigned long long x;
  744.  
  745.   x = st0 << n;
  746.  
  747.   /* Do the required multiplication and subtraction in the one operation */
  748.   asm volatile ("movl %2,%%eax; mull %4; subl %%eax,%0; sbbl %%edx,%1;
  749.                  movl %3,%%eax; mull %4; subl %%eax,%1;
  750.                  movl %2,%%eax; mull %5; subl %%eax,%1;"
  751.         :"=m" (x), "=m" (((unsigned *)&x)[1])
  752.         :"m" (st1),"m" (((unsigned *)&st1)[1]),
  753.          "m" (q),"m" (((unsigned *)&q)[1])
  754.         :"%ax","%dx");
  755.  
  756.   *y = x;
  757. }
  758.  
  759.  
  760. /* Remainder of st(0) / st(1) */
  761. /* This routine produces exact results, i.e. there is never any
  762.    rounding or truncation, etc of the result. */
  763. static void do_fprem(int round)
  764. {
  765.   FPU_REG *st1_ptr = &st(1);
  766.   char st1_tag = st1_ptr->tag;
  767.   char sign = FPU_st0_ptr->sign;
  768.  
  769.   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
  770.     {
  771.       FPU_REG tmp;
  772.       int old_cw = control_word;
  773.       int expdif = FPU_st0_ptr->exp - st1_ptr->exp;
  774.       long long q;
  775.       unsigned short saved_status;
  776.       int cc = 0;
  777.  
  778. #ifdef DENORM_OPERAND
  779.       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
  780.         (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
  781.     return;
  782. #endif DENORM_OPERAND
  783.       
  784.       /* We want the status following the denorm tests, but don't want
  785.      the status changed by the arithmetic operations. */
  786.       saved_status = partial_status;
  787.       control_word &= ~CW_RC;
  788.       control_word |= RC_CHOP;
  789.  
  790.       if (expdif < 64)
  791.     {
  792.       /* This should be the most common case */
  793.  
  794.       if ( expdif > -2 )
  795.         {
  796.           reg_div(FPU_st0_ptr, st1_ptr, &tmp, PR_64_BITS | RC_CHOP | 0x3f);
  797.  
  798.           if ( tmp.exp >= EXP_BIAS )
  799.         {
  800.           round_to_int(&tmp);  /* Fortunately, this can't overflow
  801.                       to 2^64 */
  802.           q = significand(&tmp);
  803.  
  804.           rem_kernel(significand(FPU_st0_ptr),
  805.                  &significand(&tmp),
  806.                  significand(st1_ptr),
  807.                  q, expdif);
  808.  
  809.           tmp.exp = st1_ptr->exp;
  810.         }
  811.           else
  812.         {
  813.           reg_move(FPU_st0_ptr, &tmp);
  814.           q = 0;
  815.         }
  816.           tmp.sign = sign;
  817.  
  818.           if ( (round == RC_RND) && (tmp.sigh & 0xc0000000) )
  819.         {
  820.           /* We may need to subtract st(1) once more,
  821.              to get a result <= 1/2 of st(1). */
  822.           unsigned long long x;
  823.           expdif = st1_ptr->exp - tmp.exp;
  824.           if ( expdif <= 1 )
  825.             {
  826.               if ( expdif == 0 )
  827.             x = significand(st1_ptr) - significand(&tmp);
  828.               else /* expdif is 1 */
  829.             x = (significand(st1_ptr) << 1) - significand(&tmp);
  830.               if ( (x < significand(&tmp)) ||
  831.               /* or equi-distant (from 0 & st(1)) and q is odd */
  832.               ((x == significand(&tmp)) && (q & 1) ) )
  833.             {
  834.               tmp.sign ^= (SIGN_POS^SIGN_NEG);
  835.               significand(&tmp) = x;
  836.               q++;
  837.             }
  838.             }
  839.         }
  840.  
  841.           if (q & 4) cc |= SW_C0;
  842.           if (q & 2) cc |= SW_C3;
  843.           if (q & 1) cc |= SW_C1;
  844.         }
  845.       else
  846.         {
  847.           control_word = old_cw;
  848.           setcc(0);
  849.           return;
  850.         }
  851.     }
  852.       else
  853.     {
  854.       /* There is a large exponent difference ( >= 64 ) */
  855.       /* To make much sense, the code in this section should
  856.          be done at high precision. */
  857.       int exp_1;
  858.  
  859.       /* prevent overflow here */
  860.       /* N is 'a number between 32 and 63' (p26-113) */
  861.       reg_move(FPU_st0_ptr, &tmp);
  862.       tmp.exp = EXP_BIAS + 56;
  863.       exp_1 = st1_ptr->exp;      st1_ptr->exp = EXP_BIAS;
  864.       expdif -= 56;
  865.  
  866.       reg_div(&tmp, st1_ptr, &tmp, PR_64_BITS | RC_CHOP | 0x3f);
  867.       st1_ptr->exp = exp_1;
  868.  
  869.       round_to_int(&tmp);  /* Fortunately, this can't overflow to 2^64 */
  870.  
  871.       rem_kernel(significand(FPU_st0_ptr),
  872.              &significand(&tmp),
  873.              significand(st1_ptr),
  874.              significand(&tmp),
  875.              tmp.exp - EXP_BIAS
  876.              ); 
  877.       tmp.exp = exp_1 + expdif;
  878.       tmp.sign = sign;
  879.  
  880.       /* It is possible for the operation to be complete here.
  881.          What does the IEEE standard say? The Intel 80486 manual
  882.          implies that the operation will never be completed at this
  883.          point, and the behaviour of a real 80486 confirms this.
  884.        */
  885.       if ( !(tmp.sigh | tmp.sigl) )
  886.         {
  887.           /* The result is zero */
  888.           control_word = old_cw;
  889.           partial_status = saved_status;
  890.           reg_move(&CONST_Z, FPU_st0_ptr);
  891.           FPU_st0_ptr->sign = sign;
  892. #ifdef PECULIAR_486
  893.           setcc(SW_C2);
  894. #else
  895.           setcc(0);
  896. #endif PECULIAR_486
  897.           return;
  898.         }
  899.       cc = SW_C2;
  900.     }
  901.  
  902.       control_word = old_cw;
  903.       partial_status = saved_status;
  904.       normalize_nuo(&tmp);
  905.       reg_move(&tmp, FPU_st0_ptr);
  906.       setcc(cc);
  907.  
  908.       /* The only condition to be looked for is underflow,
  909.      and it can occur here only if underflow is unmasked. */
  910.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (FPU_st0_ptr->tag != TW_Zero)
  911.       && !(control_word & CW_Underflow) )
  912.     arith_underflow(FPU_st0_ptr);
  913.  
  914.       return;
  915.     }
  916.   else if ( (FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty) )
  917.     {
  918.       stack_underflow();
  919.       return;
  920.     }
  921.   else if ( FPU_st0_tag == TW_Zero )
  922.     {
  923.       if ( st1_tag == TW_Valid )
  924.     {
  925. #ifdef DENORM_OPERAND
  926.       if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  927.         return;
  928. #endif DENORM_OPERAND
  929.  
  930.       setcc(0); return;
  931.     }
  932.       else if ( st1_tag == TW_Zero )
  933.     { arith_invalid(FPU_st0_ptr); return; } /* fprem(?,0) always invalid */
  934.       else if ( st1_tag == TW_Infinity )
  935.     { setcc(0); return; }
  936.     }
  937.   else if ( FPU_st0_tag == TW_Valid )
  938.     {
  939.       if ( st1_tag == TW_Zero )
  940.     {
  941.       arith_invalid(FPU_st0_ptr); /* fprem(Valid,Zero) is invalid */
  942.       return;
  943.     }
  944.       else if ( st1_tag != TW_NaN )
  945.     {
  946. #ifdef DENORM_OPERAND
  947.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  948.         return;
  949. #endif DENORM_OPERAND
  950.  
  951.       if ( st1_tag == TW_Infinity )
  952.         {
  953.           /* fprem(Valid,Infinity) is o.k. */
  954.           setcc(0); return;
  955.         }
  956.     }
  957.     }
  958.   else if ( FPU_st0_tag == TW_Infinity )
  959.     {
  960.       if ( st1_tag != TW_NaN )
  961.     {
  962.       arith_invalid(FPU_st0_ptr); /* fprem(Infinity,?) is invalid */
  963.       return;
  964.     }
  965.     }
  966.  
  967.   /* One of the registers must contain a NaN is we got here. */
  968.  
  969. #ifdef PARANOID
  970.   if ( (FPU_st0_tag != TW_NaN) && (st1_tag != TW_NaN) )
  971.       EXCEPTION(EX_INTERNAL | 0x118);
  972. #endif PARANOID
  973.  
  974.   real_2op_NaN(st1_ptr, FPU_st0_ptr, FPU_st0_ptr);
  975.  
  976. }
  977.  
  978.  
  979. /* ST(1) <- ST(1) * log ST;  pop ST */
  980. static void fyl2x(void)
  981. {
  982.   FPU_REG *st1_ptr = &st(1);
  983.   char st1_tag = st1_ptr->tag;
  984.  
  985.   clear_C1();
  986.   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
  987.     {
  988.       if ( FPU_st0_ptr->sign == SIGN_POS )
  989.     {
  990.       int saved_control, saved_status;
  991.  
  992. #ifdef DENORM_OPERAND
  993.       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
  994.         (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
  995.         return;
  996. #endif DENORM_OPERAND
  997.  
  998.       /* We use the general purpose arithmetic,
  999.          so we need to save these. */
  1000.       saved_status = partial_status;
  1001.       saved_control = control_word;
  1002.       control_word = FULL_PRECISION;
  1003.  
  1004.       poly_l2(FPU_st0_ptr, FPU_st0_ptr);
  1005.  
  1006.       /* Enough of the basic arithmetic is done now */
  1007.       control_word = saved_control;
  1008.       partial_status = saved_status;
  1009.  
  1010.       /* Let the multiply set the flags */
  1011.       reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
  1012.  
  1013.       pop(); FPU_st0_ptr = &st(0);
  1014.     }
  1015.       else
  1016.     {
  1017.       /* negative */
  1018.       if ( !arith_invalid(st1_ptr) )
  1019.         pop();
  1020.       return;
  1021.     }
  1022.     }
  1023.   else if ( (FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty) )
  1024.     {
  1025.       stack_underflow_pop(1);
  1026.       return;
  1027.     }
  1028.   else if ( (FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN) )
  1029.     {
  1030.       if ( !real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr) )
  1031.     pop();
  1032.       return;
  1033.     }
  1034.   else if ( (FPU_st0_tag <= TW_Zero) && (st1_tag <= TW_Zero) )
  1035.     {
  1036.       /* one of the args is zero, the other valid, or both zero */
  1037.       if ( FPU_st0_tag == TW_Zero )
  1038.     {
  1039.       if ( st1_tag == TW_Zero )
  1040.         {
  1041.           /* Both args zero is invalid */
  1042.           if ( !arith_invalid(st1_ptr) )
  1043.         pop();
  1044.         }
  1045. #ifdef PECULIAR_486
  1046.       /* This case is not specifically covered in the manual,
  1047.          but divide-by-zero would seem to be the best response.
  1048.          However, a real 80486 does it this way... */
  1049.       else if ( FPU_st0_ptr->tag == TW_Infinity )
  1050.         {
  1051.           reg_move(&CONST_INF, st1_ptr);
  1052.           pop();
  1053.         }
  1054. #endif PECULIAR_486
  1055.       else
  1056.         {
  1057.           if ( !divide_by_zero(st1_ptr->sign^SIGN_NEG^SIGN_POS, st1_ptr) )
  1058.         pop();
  1059.         }
  1060.       return;
  1061.     }
  1062.       else
  1063.     {
  1064.       /* st(1) contains zero, st(0) valid <> 0 */
  1065.       /* Zero is the valid answer */
  1066.       char sign = st1_ptr->sign;
  1067.  
  1068.       if ( FPU_st0_ptr->sign == SIGN_NEG )
  1069.         {
  1070.           /* log(negative) */
  1071.           if ( !arith_invalid(st1_ptr) )
  1072.         pop();
  1073.           return;
  1074.         }
  1075.  
  1076. #ifdef DENORM_OPERAND
  1077.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1078.         return;
  1079. #endif DENORM_OPERAND
  1080.  
  1081.       if ( FPU_st0_ptr->exp < EXP_BIAS ) sign ^= SIGN_NEG^SIGN_POS;
  1082.       pop(); FPU_st0_ptr = &st(0);
  1083.       reg_move(&CONST_Z, FPU_st0_ptr);
  1084.       FPU_st0_ptr->sign = sign;
  1085.       return;
  1086.     }
  1087.     }
  1088.   /* One or both arg must be an infinity */
  1089.   else if ( FPU_st0_tag == TW_Infinity )
  1090.     {
  1091.       if ( (FPU_st0_ptr->sign == SIGN_NEG) || (st1_tag == TW_Zero) )
  1092.     {
  1093.       /* log(-infinity) or 0*log(infinity) */
  1094.       if ( !arith_invalid(st1_ptr) )
  1095.         pop();
  1096.       return;
  1097.     }
  1098.       else
  1099.     {
  1100.       char sign = st1_ptr->sign;
  1101.  
  1102. #ifdef DENORM_OPERAND
  1103.       if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1104.         return;
  1105. #endif DENORM_OPERAND
  1106.  
  1107.       pop(); FPU_st0_ptr = &st(0);
  1108.       reg_move(&CONST_INF, FPU_st0_ptr);
  1109.       FPU_st0_ptr->sign = sign;
  1110.       return;
  1111.     }
  1112.     }
  1113.   /* st(1) must be infinity here */
  1114.   else if ( (FPU_st0_tag == TW_Valid) && (FPU_st0_ptr->sign == SIGN_POS) )
  1115.     {
  1116.       if ( FPU_st0_ptr->exp >= EXP_BIAS )
  1117.     {
  1118.       if ( (FPU_st0_ptr->exp == EXP_BIAS) &&
  1119.           (FPU_st0_ptr->sigh == 0x80000000) &&
  1120.           (FPU_st0_ptr->sigl == 0) )
  1121.         {
  1122.           /* st(0) holds 1.0 */
  1123.           /* infinity*log(1) */
  1124.           if ( !arith_invalid(st1_ptr) )
  1125.         pop();
  1126.           return;
  1127.         }
  1128.       /* st(0) is positive and > 1.0 */
  1129.       pop();
  1130.     }
  1131.       else
  1132.     {
  1133.       /* st(0) is positive and < 1.0 */
  1134.  
  1135. #ifdef DENORM_OPERAND
  1136.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1137.         return;
  1138. #endif DENORM_OPERAND
  1139.  
  1140.       st1_ptr->sign ^= SIGN_NEG;
  1141.       pop();
  1142.     }
  1143.       return;
  1144.     }
  1145.   else
  1146.     {
  1147.       /* st(0) must be zero or negative */
  1148.       if ( FPU_st0_ptr->tag == TW_Zero )
  1149.     {
  1150.       /* This should be invalid, but a real 80486 is happy with it. */
  1151. #ifndef PECULIAR_486
  1152.       if ( !divide_by_zero(st1_ptr->sign, st1_ptr) )
  1153. #endif PECULIAR_486
  1154.         {
  1155.           st1_ptr->sign ^= SIGN_NEG^SIGN_POS;
  1156.           pop();
  1157.         }
  1158.     }
  1159.       else
  1160.     {
  1161.       /* log(negative) */
  1162.       if ( !arith_invalid(st1_ptr) )
  1163.         pop();
  1164.     }
  1165.       return;
  1166.     }
  1167. }
  1168.  
  1169.  
  1170. static void fpatan(void)
  1171. {
  1172.   FPU_REG *st1_ptr = &st(1);
  1173.   char st1_tag = st1_ptr->tag;
  1174.   char st1_sign = st1_ptr->sign, st0_sign = FPU_st0_ptr->sign;
  1175.  
  1176.   clear_C1();
  1177.   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
  1178.     {
  1179.       int saved_control, saved_status;
  1180.       FPU_REG sum;
  1181.       char inverted;
  1182.  
  1183. #ifdef DENORM_OPERAND
  1184.       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
  1185.         (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
  1186.     return;
  1187. #endif DENORM_OPERAND
  1188.  
  1189.       /* We use the general purpose arithmetic so we need to save these. */
  1190.       saved_status = partial_status;
  1191.       saved_control = control_word;
  1192.       control_word = FULL_PRECISION;
  1193.  
  1194.       st1_ptr->sign = FPU_st0_ptr->sign = SIGN_POS;
  1195.       if ( (compare(st1_ptr) & ~COMP_Denormal) == COMP_A_lt_B )
  1196.     {
  1197.       inverted = 1;
  1198.       reg_div(FPU_st0_ptr, st1_ptr, &sum, FULL_PRECISION);
  1199.     }
  1200.       else
  1201.     {
  1202.       inverted = 0;
  1203.       if ( (st0_sign == 0) &&
  1204.           (st1_ptr->exp - FPU_st0_ptr->exp < -64) )
  1205.         {
  1206.           control_word = saved_control;
  1207.           partial_status = saved_status;
  1208.           reg_div(st1_ptr, FPU_st0_ptr, st1_ptr,
  1209.               control_word | PR_64_BITS);
  1210.           st1_ptr->sign = st1_sign;
  1211.           pop();
  1212.           set_precision_flag_down();
  1213.           return;
  1214.         }
  1215.       reg_div(st1_ptr, FPU_st0_ptr, &sum, FULL_PRECISION);
  1216.     }
  1217.  
  1218.       poly_atan(&sum);
  1219.  
  1220.       if ( inverted )
  1221.     {
  1222.       reg_sub(&CONST_PI2, &sum, &sum, FULL_PRECISION);
  1223.     }
  1224.       if ( st0_sign )
  1225.     {
  1226.       reg_sub(&CONST_PI, &sum, &sum, FULL_PRECISION);
  1227.     }
  1228.       sum.sign = st1_sign;
  1229.  
  1230.       /* All of the basic arithmetic is done now */
  1231.       control_word = saved_control;
  1232.       partial_status = saved_status;
  1233.  
  1234.       reg_move(&sum, st1_ptr);
  1235.     }
  1236.   else if ( (FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty) )
  1237.     {
  1238.       stack_underflow_pop(1);
  1239.       return;
  1240.     }
  1241.   else if ( (FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN) )
  1242.     {
  1243.       if ( !real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr) )
  1244.       pop();
  1245.       return;
  1246.     }
  1247.   else if ( (FPU_st0_tag == TW_Infinity) || (st1_tag == TW_Infinity) )
  1248.     {
  1249.       char sign = st1_ptr->sign;
  1250.       if ( FPU_st0_tag == TW_Infinity )
  1251.     {
  1252.       if ( st1_tag == TW_Infinity )
  1253.         {
  1254.           if ( FPU_st0_ptr->sign == SIGN_POS )
  1255.         { reg_move(&CONST_PI4, st1_ptr); }
  1256.           else
  1257.         reg_add(&CONST_PI4, &CONST_PI2, st1_ptr, FULL_PRECISION);
  1258.         }
  1259.       else
  1260.         {
  1261. #ifdef DENORM_OPERAND
  1262.           if ( st1_tag != TW_Zero )
  1263.         {
  1264.           if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1265.             return;
  1266.         }
  1267. #endif DENORM_OPERAND
  1268.  
  1269.           if ( FPU_st0_ptr->sign == SIGN_POS )
  1270.         {
  1271.           reg_move(&CONST_Z, st1_ptr);
  1272.           st1_ptr->sign = sign;   /* An 80486 preserves the sign */
  1273.           pop();
  1274.           return;
  1275.         }
  1276.           else
  1277.         reg_move(&CONST_PI, st1_ptr);
  1278.         }
  1279.     }
  1280.       else
  1281.     {
  1282.       /* st(1) is infinity, st(0) not infinity */
  1283. #ifdef DENORM_OPERAND
  1284.       if ( FPU_st0_tag != TW_Zero )
  1285.         {
  1286.           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1287.         return;
  1288.         }
  1289. #endif DENORM_OPERAND
  1290.  
  1291.       reg_move(&CONST_PI2, st1_ptr);
  1292.     }
  1293.       st1_ptr->sign = sign;
  1294.     }
  1295.   else if ( st1_tag == TW_Zero )
  1296.     {
  1297.       /* st(0) must be valid or zero */
  1298.       char sign = st1_ptr->sign;
  1299.  
  1300. #ifdef DENORM_OPERAND
  1301.       if ( FPU_st0_tag != TW_Zero )
  1302.     {
  1303.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1304.         return;
  1305.     }
  1306. #endif DENORM_OPERAND
  1307.  
  1308.       if ( FPU_st0_ptr->sign == SIGN_POS )
  1309.     { /* An 80486 preserves the sign */ pop(); return; }
  1310.       else
  1311.     reg_move(&CONST_PI, st1_ptr);
  1312.       st1_ptr->sign = sign;
  1313.     }
  1314.   else if ( FPU_st0_tag == TW_Zero )
  1315.     {
  1316.       /* st(1) must be TW_Valid here */
  1317.       char sign = st1_ptr->sign;
  1318.  
  1319. #ifdef DENORM_OPERAND
  1320.       if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1321.     return;
  1322. #endif DENORM_OPERAND
  1323.  
  1324.       reg_move(&CONST_PI2, st1_ptr);
  1325.       st1_ptr->sign = sign;
  1326.     }
  1327. #ifdef PARANOID
  1328.   else
  1329.     EXCEPTION(EX_INTERNAL | 0x125);
  1330. #endif PARANOID
  1331.  
  1332.   pop();
  1333.   set_precision_flag_up();  /* We do not really know if up or down */
  1334. }
  1335.  
  1336.  
  1337. static void fprem(void)
  1338. {
  1339.   do_fprem(RC_CHOP);
  1340. }
  1341.  
  1342.  
  1343. static void fprem1(void)
  1344. {
  1345.   do_fprem(RC_RND);
  1346. }
  1347.  
  1348.  
  1349. static void fyl2xp1(void)
  1350. {
  1351.   FPU_REG *st1_ptr = &st(1);
  1352.   char st1_tag = st1_ptr->tag;
  1353.  
  1354.   clear_C1();
  1355.   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
  1356.     {
  1357.       int saved_control, saved_status;
  1358.  
  1359. #ifdef DENORM_OPERAND
  1360.       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
  1361.         (st1_ptr->exp <= EXP_UNDER)) && denormal_operand() )
  1362.     return;
  1363. #endif DENORM_OPERAND
  1364.  
  1365.       /* We use the general purpose arithmetic so we need to save these. */
  1366.       saved_status = partial_status;
  1367.       saved_control = control_word;
  1368.       control_word = FULL_PRECISION;
  1369.  
  1370.       if ( poly_l2p1(FPU_st0_ptr, FPU_st0_ptr) )
  1371.     {
  1372. #ifdef PECULIAR_486   /* Stupid 80486 doesn't worry about log(negative). */
  1373.       st1_ptr->sign ^= SIGN_POS^SIGN_NEG;
  1374.       control_word = saved_control;
  1375.       partial_status = saved_status;
  1376.       set_precision_flag_down();
  1377. #else
  1378.       if ( arith_invalid(st1_ptr) )  /* poly_l2p1() returned invalid */
  1379.         return;
  1380. #endif PECULIAR_486
  1381.       pop(); return;
  1382.     }
  1383.       
  1384.       /* Enough of the basic arithmetic is done now */
  1385.       control_word = saved_control;
  1386.       partial_status = saved_status;
  1387.  
  1388.       /* Let the multiply set the flags */
  1389.       reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
  1390.  
  1391.       pop();
  1392.     }
  1393.   else if ( (FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty) )
  1394.     {
  1395.       stack_underflow_pop(1);
  1396.       return;
  1397.     }
  1398.   else if ( FPU_st0_tag == TW_Zero )
  1399.     {
  1400.       if ( st1_tag <= TW_Zero )
  1401.     {
  1402. #ifdef DENORM_OPERAND
  1403.       if ( (st1_tag == TW_Valid) && (st1_ptr->exp <= EXP_UNDER) &&
  1404.           (denormal_operand()) )
  1405.         return;
  1406. #endif DENORM_OPERAND
  1407.       
  1408.       FPU_st0_ptr->sign ^= st1_ptr->sign;
  1409.       reg_move(FPU_st0_ptr, st1_ptr);
  1410.     }
  1411.       else if ( st1_tag == TW_Infinity )
  1412.     {
  1413.       /* Infinity*log(1) */
  1414.       if ( !arith_invalid(st1_ptr) )
  1415.         pop();
  1416.       return;
  1417.     }
  1418.       else if ( st1_tag == TW_NaN )
  1419.     {
  1420.       if ( !real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr) )
  1421.         pop();
  1422.       return;
  1423.     }
  1424. #ifdef PARANOID
  1425.       else
  1426.     {
  1427.       EXCEPTION(EX_INTERNAL | 0x116);
  1428.       return;
  1429.     }
  1430. #endif PARANOID
  1431.       pop(); return;
  1432.     }
  1433.   else if ( FPU_st0_tag == TW_Valid )
  1434.     {
  1435.       if ( st1_tag == TW_Zero )
  1436.     {
  1437.       if ( FPU_st0_ptr->sign == SIGN_NEG )
  1438.         {
  1439.           if ( FPU_st0_ptr->exp >= EXP_BIAS )
  1440.         {
  1441.           /* st(0) holds <= -1.0 */
  1442. #ifdef PECULIAR_486   /* Stupid 80486 doesn't worry about log(negative). */
  1443.           st1_ptr->sign ^= SIGN_POS^SIGN_NEG;
  1444. #else
  1445.           if ( arith_invalid(st1_ptr) ) return;
  1446. #endif PECULIAR_486
  1447.           pop(); return;
  1448.         }
  1449. #ifdef DENORM_OPERAND
  1450.           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1451.         return;
  1452. #endif DENORM_OPERAND
  1453.           st1_ptr->sign ^= SIGN_POS^SIGN_NEG;
  1454.           pop(); return;
  1455.         }
  1456. #ifdef DENORM_OPERAND
  1457.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1458.         return;
  1459. #endif DENORM_OPERAND
  1460.       pop(); return;
  1461.     }
  1462.       if ( st1_tag == TW_Infinity )
  1463.     {
  1464.       if ( FPU_st0_ptr->sign == SIGN_NEG )
  1465.         {
  1466.           if ( (FPU_st0_ptr->exp >= EXP_BIAS) &&
  1467.           !((FPU_st0_ptr->sigh == 0x80000000) &&
  1468.             (FPU_st0_ptr->sigl == 0)) )
  1469.         {
  1470.           /* st(0) holds < -1.0 */
  1471. #ifdef PECULIAR_486   /* Stupid 80486 doesn't worry about log(negative). */
  1472.           st1_ptr->sign ^= SIGN_POS^SIGN_NEG;
  1473. #else
  1474.           if ( arith_invalid(st1_ptr) ) return;
  1475. #endif PECULIAR_486
  1476.           pop(); return;
  1477.         }
  1478. #ifdef DENORM_OPERAND
  1479.           if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1480.         return;
  1481. #endif DENORM_OPERAND
  1482.           st1_ptr->sign ^= SIGN_POS^SIGN_NEG;
  1483.           pop(); return;
  1484.         }
  1485. #ifdef DENORM_OPERAND
  1486.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1487.         return;
  1488. #endif DENORM_OPERAND
  1489.       pop(); return;
  1490.     }
  1491.       if ( st1_tag == TW_NaN )
  1492.     {
  1493.       if ( !real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr) )
  1494.         pop();
  1495.       return;
  1496.     }
  1497.     }
  1498.   else if ( FPU_st0_tag == TW_NaN )
  1499.     {
  1500.       if ( !real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr) )
  1501.     pop();
  1502.       return;
  1503.     }
  1504.   else if ( FPU_st0_tag == TW_Infinity )
  1505.     {
  1506.       if ( st1_tag == TW_NaN )
  1507.     {
  1508.       if ( !real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr) )
  1509.         pop();
  1510.       return;
  1511.     }
  1512.       else if ( FPU_st0_ptr->sign == SIGN_NEG )
  1513.     {
  1514.       int exponent = st1_ptr->exp;
  1515. #ifndef PECULIAR_486
  1516.       /* This should have higher priority than denormals, but... */
  1517.       if ( arith_invalid(st1_ptr) )  /* log(-infinity) */
  1518.         return;
  1519. #endif PECULIAR_486
  1520. #ifdef DENORM_OPERAND
  1521.       if ( st1_tag != TW_Zero )
  1522.         {
  1523.           if ( (exponent <= EXP_UNDER) && (denormal_operand()) )
  1524.         return;
  1525.         }
  1526. #endif DENORM_OPERAND
  1527. #ifdef PECULIAR_486
  1528.       /* Denormal operands actually get higher priority */
  1529.       if ( arith_invalid(st1_ptr) )  /* log(-infinity) */
  1530.         return;
  1531. #endif PECULIAR_486
  1532.       pop();
  1533.       return;
  1534.     }
  1535.       else if ( st1_tag == TW_Zero )
  1536.     {
  1537.       /* log(infinity) */
  1538.       if ( !arith_invalid(st1_ptr) )
  1539.         pop();
  1540.       return;
  1541.     }
  1542.     
  1543.       /* st(1) must be valid here. */
  1544.  
  1545. #ifdef DENORM_OPERAND
  1546.       if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1547.     return;
  1548. #endif DENORM_OPERAND
  1549.  
  1550.       /* The Manual says that log(Infinity) is invalid, but a real
  1551.      80486 sensibly says that it is o.k. */
  1552.       { char sign = st1_ptr->sign;
  1553.     reg_move(&CONST_INF, st1_ptr);
  1554.     st1_ptr->sign = sign;
  1555.       }
  1556.       pop();
  1557.       return;
  1558.     }
  1559. #ifdef PARANOID
  1560.   else
  1561.     {
  1562.       EXCEPTION(EX_INTERNAL | 0x117);
  1563.     }
  1564. #endif PARANOID
  1565. }
  1566.  
  1567.  
  1568. static void fscale(void)
  1569. {
  1570.   FPU_REG *st1_ptr = &st(1);
  1571.   char st1_tag = st1_ptr->tag;
  1572.   int old_cw = control_word;
  1573.   char sign = FPU_st0_ptr->sign;
  1574.  
  1575.   clear_C1();
  1576.   if ( !((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid)) )
  1577.     {
  1578.       long scale;
  1579.       FPU_REG tmp;
  1580.  
  1581. #ifdef DENORM_OPERAND
  1582.       if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
  1583.         (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
  1584.     return;
  1585. #endif DENORM_OPERAND
  1586.  
  1587.       if ( st1_ptr->exp > EXP_BIAS + 30 )
  1588.     {
  1589.       /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
  1590.       char sign;
  1591.  
  1592.       if ( st1_ptr->sign == SIGN_POS )
  1593.         {
  1594.           EXCEPTION(EX_Overflow);
  1595.           sign = FPU_st0_ptr->sign;
  1596.           reg_move(&CONST_INF, FPU_st0_ptr);
  1597.           FPU_st0_ptr->sign = sign;
  1598.         }
  1599.       else
  1600.         {
  1601.           EXCEPTION(EX_Underflow);
  1602.           sign = FPU_st0_ptr->sign;
  1603.           reg_move(&CONST_Z, FPU_st0_ptr);
  1604.           FPU_st0_ptr->sign = sign;
  1605.         }
  1606.       return;
  1607.     }
  1608.  
  1609.       control_word &= ~CW_RC;
  1610.       control_word |= RC_CHOP;
  1611.       reg_move(st1_ptr, &tmp);
  1612.       round_to_int(&tmp);               /* This can never overflow here */
  1613.       control_word = old_cw;
  1614.       scale = st1_ptr->sign ? -tmp.sigl : tmp.sigl;
  1615.       scale += FPU_st0_ptr->exp;
  1616.       FPU_st0_ptr->exp = scale;
  1617.  
  1618.       /* Use round_reg() to properly detect under/overflow etc */
  1619.       round_reg(FPU_st0_ptr, 0, control_word);
  1620.  
  1621.       return;
  1622.     }
  1623.   else if ( FPU_st0_tag == TW_Valid )
  1624.     {
  1625.       if ( st1_tag == TW_Zero )
  1626.     {
  1627.  
  1628. #ifdef DENORM_OPERAND
  1629.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1630.         return;
  1631. #endif DENORM_OPERAND
  1632.  
  1633.       return;
  1634.     }
  1635.       if ( st1_tag == TW_Infinity )
  1636.     {
  1637. #ifdef DENORM_OPERAND
  1638.       if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1639.         return;
  1640. #endif DENORM_OPERAND
  1641.  
  1642.       if ( st1_ptr->sign == SIGN_POS )
  1643.         { reg_move(&CONST_INF, FPU_st0_ptr); }
  1644.       else
  1645.           reg_move(&CONST_Z, FPU_st0_ptr);
  1646.       FPU_st0_ptr->sign = sign;
  1647.       return;
  1648.     }
  1649.       if ( st1_tag == TW_NaN )
  1650.     { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
  1651.     }
  1652.   else if ( FPU_st0_tag == TW_Zero )
  1653.     {
  1654.       if ( st1_tag == TW_Valid )
  1655.     {
  1656.  
  1657. #ifdef DENORM_OPERAND
  1658.       if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1659.         return;
  1660. #endif DENORM_OPERAND
  1661.  
  1662.       return;
  1663.     }
  1664.       else if ( st1_tag == TW_Zero ) { return; }
  1665.       else if ( st1_tag == TW_Infinity )
  1666.     {
  1667.       if ( st1_ptr->sign == SIGN_NEG )
  1668.         return;
  1669.       else
  1670.         {
  1671.           arith_invalid(FPU_st0_ptr); /* Zero scaled by +Infinity */
  1672.           return;
  1673.         }
  1674.     }
  1675.       else if ( st1_tag == TW_NaN )
  1676.     { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
  1677.     }
  1678.   else if ( FPU_st0_tag == TW_Infinity )
  1679.     {
  1680.       if ( st1_tag == TW_Valid )
  1681.     {
  1682.  
  1683. #ifdef DENORM_OPERAND
  1684.       if ( (st1_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
  1685.         return;
  1686. #endif DENORM_OPERAND
  1687.  
  1688.       return;
  1689.     }
  1690.       if ( ((st1_tag == TW_Infinity) && (st1_ptr->sign == SIGN_POS))
  1691.       || (st1_tag == TW_Zero) )
  1692.     return;
  1693.       else if ( st1_tag == TW_Infinity )
  1694.     {
  1695.       arith_invalid(FPU_st0_ptr); /* Infinity scaled by -Infinity */
  1696.       return;
  1697.     }
  1698.       else if ( st1_tag == TW_NaN )
  1699.     { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
  1700.     }
  1701.   else if ( FPU_st0_tag == TW_NaN )
  1702.     {
  1703.       if ( st1_tag != TW_Empty )
  1704.     { real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr); return; }
  1705.     }
  1706.  
  1707. #ifdef PARANOID
  1708.   if ( !((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty)) )
  1709.     {
  1710.       EXCEPTION(EX_INTERNAL | 0x115);
  1711.       return;
  1712.     }
  1713. #endif
  1714.  
  1715.   /* At least one of st(0), st(1) must be empty */
  1716.   stack_underflow();
  1717.  
  1718. }
  1719.  
  1720.  
  1721. /*---------------------------------------------------------------------------*/
  1722.  
  1723. static FUNC const trig_table_a[] = {
  1724.   f2xm1, fyl2x, fptan, fpatan, fxtract, fprem1, fdecstp, fincstp
  1725. };
  1726.  
  1727. void trig_a(void)
  1728. {
  1729.   (trig_table_a[FPU_rm])();
  1730. }
  1731.  
  1732.  
  1733. static FUNC const trig_table_b[] =
  1734.   {
  1735.     fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, fsin, fcos
  1736.   };
  1737.  
  1738. void trig_b(void)
  1739. {
  1740.   (trig_table_b[FPU_rm])();
  1741. }
  1742.